Calculating the influence from the nearest feature in both R and GRASS GIS using the oneimpact pacakge

# libraries
library(rgrass7)
library(raster)
library(terra)

library(NinaR)
library(dplyr)

library(viridis)

library(oneimpact)

Connect to GRASS.

# start GRASS and connect to my mapset
grassdir <- system("grass78 --config path", intern = T)
gisDB <- "/data/grass"
loc <- "ETRS_33N/"
ms <- "u_bb_cuminf"
# initGRASS(gisBase = grassdir,
#           home = tempdir(), 
#           override = T,
#           gisDbase = gisDB,
#           location = loc, 
#           mapset = ms)

# more directly within NINA
NinaR::grassConnect(mapset = ms)

## gisdbase    /data/grass 
## location    ETRS_33N 
## mapset      u_bb_cuminf 
## rows        361 
## columns     478 
## north       6658900 
## south       6622800 
## west        146900 
## east        194700 
## nsres       100 
## ewres       100

# copy test region vector and map of private cabins

# check region
gmeta() # g.copy is not affected by the region

## gisdbase    /data/grass 
## location    ETRS_33N 
## mapset      u_bb_cuminf 
## rows        361 
## columns     478 
## north       6658900 
## south       6622800 
## west        146900 
## east        194700 
## nsres       100 
## ewres       100

# copy region vector
rgrass7::execGRASS("g.copy", parameters = list(vector = "region_test_influence@u_bernardo.brandao,region_test_influence"),
                   flags = c("overwrite"))

## Warning in rgrass7::execGRASS("g.copy", parameters = list(vector = "region_test_influence@u_bernardo.brandao,region_test_influence"), : The command:
## g.copy --overwrite vector=region_test_influence@u_bernardo.brandao,region_test_influence
## produced at least one warning during execution:
## Copy vector <region_test_influence@u_bernardo.brandao> to current mapset as
## <region_test_influence>
## WARNING: Vector map <region_test_influence> already exists and will be
##          overwritten

## Copy vector <region_test_influence@u_bernardo.brandao> to current mapset as
## <region_test_influence>
## WARNING: Vector map <region_test_influence> already exists and will be
##          overwritten

# copy private cabins
rgrass7::execGRASS("g.copy", parameters = list(raster = "private_cabins_rast@u_bernardo.brandao,private_cabins_rast"),
                   flags = c("overwrite"))

## Copy raster <private_cabins_rast@u_bernardo.brandao> to current mapset as
## <private_cabins_rast>

# define region

# I defined a region_test_influence through the GRASS GUI and using v.in.region output=region_test_influence
rgrass7::execGRASS("g.region", parameters = list(vector = "region_test_influence", align = "private_cabins_rast"),
          flags = "p")

## projection: 1 (UTM)
## zone:       33
## datum:      etrs89
## ellipsoid:  grs80
## north:      6658900
## south:      6622800
## west:       146900
## east:       194700
## nsres:      100
## ewres:      100
## rows:       361
## cols:       478
## cells:      172558

# create a new small map for the current region, just to ease visualization in R
rgrass7::execGRASS("r.mapcalc", expression = "private_cabins_sub = if(!isnull(private_cabins_rast), 1, null())",
          flags = c("overwrite", "quiet"))

# show input map and region here
rgrass7::use_sp()
cabins <- readRAST(c("private_cabins_rast_sub")) %>% 
  raster::raster() %>% 
  terra::rast()

## Warning in getProjectionRef(x, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = OVERRIDE_PROJ_DATUM_WITH_TOWGS84, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
##  but +towgs84= values preserved

## Warning in rgdal::GDALinfo(fname, silent = ignore.stderr): statistics not
## supported by this driver

## /data/grass/ETRS_33N/u_bb_cuminf/cellhd/private_cabins_rast_sub has GDAL driver GRASS 
## and has 361 rows and 478 columns

## Warning in getProjectionRef(x, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = OVERRIDE_PROJ_DATUM_WITH_TOWGS84, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
##  but +towgs84= values preserved

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on GRS80 ellipsoid in Proj4
## definition

terra::plot(cabins, col = "black")

# Euclidean distance

# R
cabins_dist_r <- calc_influence_nearest(cabins, where = "R")

# GRASS
name_var <- "private_cabins_sub"
cabins_dist_names <- calc_influence_nearest(name_var, where = "GRASS", quiet = F, overwrite = T)

## projection: 1 (UTM)
## zone:       33
## datum:      etrs89
## ellipsoid:  grs80
## north:      6658900
## south:      6622800
## west:       146900
## east:       194700
## nsres:      100
## ewres:      100
## rows:       361
## cols:       478
## cells:      172558
## Calculating Euclidean distance...
## Reading raster map <private_cabins_sub>...
##    0%   3%   6%   9%  12%  15%  18%  21%  24%  27%  30%  33%  36%  39%  42%  45%  48%  51%  54%  57%  60%  63%  66%  69%  72%  75%  78%  81%  84%  87%  90%  93%  96%  99% 100%
## Writing output raster maps...
##    0%   3%   6%   9%  12%  15%  18%  21%  24%  27%  30%  33%  36%  39%  42%  45%  48%  51%  54%  57%  60%  63%  66%  69%  72%  75%  78%  81%  84%  87%  90%  93%  96%  99% 100%

# visualize
cabins_dist_grass <- readRAST(cabins_dist_names) %>% 
  raster::raster() %>% 
  terra::rast()

## Warning in getProjectionRef(x, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = OVERRIDE_PROJ_DATUM_WITH_TOWGS84, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
##  but +towgs84= values preserved

## Warning in rgdal::GDALinfo(fname, silent = ignore.stderr): statistics not
## supported by this driver

## /data/grass/ETRS_33N/u_bb_cuminf/cellhd/private_cabins_sub_inf_nearest has GDAL driver GRASS 
## and has 361 rows and 478 columns

## Warning in getProjectionRef(x, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = OVERRIDE_PROJ_DATUM_WITH_TOWGS84, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
##  but +towgs84= values preserved

terra::plot(c(cabins_dist_r, cabins_dist_grass), 
            main = paste("Euclidean distance from private cabins - ", 
                         c("R", "GRASS")))

# Log distance

# R
cabins_logdist_r <- calc_influence_nearest(cabins, transform = "log", log_base = 10, where = "R")

# GRASS
name_var <- "private_cabins_sub"
cabins_logdist_names <- calc_influence_nearest(name_var, transform = "log", log_base = 10, 
                                               where = "GRASS", quiet = F, overwrite = T)

## projection: 1 (UTM)
## zone:       33
## datum:      etrs89
## ellipsoid:  grs80
## north:      6658900
## south:      6622800
## west:       146900
## east:       194700
## nsres:      100
## ewres:      100
## rows:       361
## cols:       478
## cells:      172558
## Calculating Euclidean distance...
## Reading raster map <private_cabins_sub>...
##    0%   3%   6%   9%  12%  15%  18%  21%  24%  27%  30%  33%  36%  39%  42%  45%  48%  51%  54%  57%  60%  63%  66%  69%  72%  75%  78%  81%  84%  87%  90%  93%  96%  99% 100%
## Writing output raster maps...
##    0%   3%   6%   9%  12%  15%  18%  21%  24%  27%  30%  33%  36%  39%  42%  45%  48%  51%  54%  57%  60%  63%  66%  69%  72%  75%  78%  81%  84%  87%  90%  93%  96%  99% 100%
## [1] "log(A + 1.000000, 10.000000)"
## Calculating log-distance...
## projection: 1 (UTM)
## zone:       33
## datum:      etrs89
## ellipsoid:  grs80
## north:      6658900
## south:      6622800
## west:       146900
## east:       194700
## nsres:      100
## ewres:      100
## rows:       361
## cols:       478
## cells:      172558

# visualize
cabins_logdist_grass <- readRAST(cabins_logdist_names) %>% 
  raster::raster() %>% 
  terra::rast()

## Warning in getProjectionRef(x, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = OVERRIDE_PROJ_DATUM_WITH_TOWGS84, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
##  but +towgs84= values preserved

## Warning in rgdal::GDALinfo(fname, silent = ignore.stderr): statistics not
## supported by this driver

## /data/grass/ETRS_33N/u_bb_cuminf/cellhd/private_cabins_sub_inf_nearest_log has GDAL driver GRASS 
## and has 361 rows and 478 columns

## Warning in getProjectionRef(x, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = OVERRIDE_PROJ_DATUM_WITH_TOWGS84, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
##  but +towgs84= values preserved

terra::plot(c(cabins_logdist_r, cabins_logdist_grass), 
            main = paste("Log-distance from private cabins - ", 
                         c("R", "GRASS")))

Square root distance

# R
cabins_sqrtdist_r <- calc_influence_nearest(cabins, transform = "sqrt", where = "R")

# GRASS
name_var <- "private_cabins_sub"
cabins_sqrtdist_names <- calc_influence_nearest(name_var, transform = "sqrt", 
                                               where = "GRASS", quiet = F, overwrite = T)

## projection: 1 (UTM)
## zone:       33
## datum:      etrs89
## ellipsoid:  grs80
## north:      6658900
## south:      6622800
## west:       146900
## east:       194700
## nsres:      100
## ewres:      100
## rows:       361
## cols:       478
## cells:      172558
## Calculating Euclidean distance...
## Reading raster map <private_cabins_sub>...
##    0%   3%   6%   9%  12%  15%  18%  21%  24%  27%  30%  33%  36%  39%  42%  45%  48%  51%  54%  57%  60%  63%  66%  69%  72%  75%  78%  81%  84%  87%  90%  93%  96%  99% 100%
## Writing output raster maps...
##    0%   3%   6%   9%  12%  15%  18%  21%  24%  27%  30%  33%  36%  39%  42%  45%  48%  51%  54%  57%  60%  63%  66%  69%  72%  75%  78%  81%  84%  87%  90%  93%  96%  99% 100%
## [1] "sqrt(A + 1.000000)"
## Calculating sqrt-distance...
## projection: 1 (UTM)
## zone:       33
## datum:      etrs89
## ellipsoid:  grs80
## north:      6658900
## south:      6622800
## west:       146900
## east:       194700
## nsres:      100
## ewres:      100
## rows:       361
## cols:       478
## cells:      172558

# visualize
cabins_sqrtdist_grass <- readRAST(cabins_sqrtdist_names) %>% 
  raster::raster() %>% 
  terra::rast()

## Warning in getProjectionRef(x, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = OVERRIDE_PROJ_DATUM_WITH_TOWGS84, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
##  but +towgs84= values preserved

## Warning in rgdal::GDALinfo(fname, silent = ignore.stderr): statistics not
## supported by this driver

## /data/grass/ETRS_33N/u_bb_cuminf/cellhd/private_cabins_sub_inf_nearest_sqrt has GDAL driver GRASS 
## and has 361 rows and 478 columns

## Warning in getProjectionRef(x, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = OVERRIDE_PROJ_DATUM_WITH_TOWGS84, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
##  but +towgs84= values preserved

terra::plot(c(cabins_sqrtdist_r, cabins_sqrtdist_grass), 
            main = paste("Sqrt-distance from private cabins - ", 
                         c("R", "GRASS")))

Exponental decay influence

Example with ZoI = 1000m

# R
cabins_exp_decay_r <- calc_influence_nearest(cabins, transform = "exp_decay", 
                                             zoi = 1000, where = "R")

# GRASS
name_var <- "private_cabins_sub"
cabins_exp_decay_names <- calc_influence_nearest(name_var, transform = "exp_decay",
                                                 zoi = 1000, where = "GRASS", 
                                                 quiet = F, overwrite = T)

## projection: 1 (UTM)
## zone:       33
## datum:      etrs89
## ellipsoid:  grs80
## north:      6658900
## south:      6622800
## west:       146900
## east:       194700
## nsres:      100
## ewres:      100
## rows:       361
## cols:       478
## cells:      172558
## Calculating Euclidean distance...
## Reading raster map <private_cabins_sub>...
##    0%   3%   6%   9%  12%  15%  18%  21%  24%  27%  30%  33%  36%  39%  42%  45%  48%  51%  54%  57%  60%  63%  66%  69%  72%  75%  78%  81%  84%  87%  90%  93%  96%  99% 100%
## Writing output raster maps...
##    0%   3%   6%   9%  12%  15%  18%  21%  24%  27%  30%  33%  36%  39%  42%  45%  48%  51%  54%  57%  60%  63%  66%  69%  72%  75%  78%  81%  84%  87%  90%  93%  96%  99% 100%
## [1] "1.000000 * exp(-0.002773 * A)"
## Calculating exponential decay influence...
## projection: 1 (UTM)
## zone:       33
## datum:      etrs89
## ellipsoid:  grs80
## north:      6658900
## south:      6622800
## west:       146900
## east:       194700
## nsres:      100
## ewres:      100
## rows:       361
## cols:       478
## cells:      172558

# visualize
cabins_exp_decay_grass <- readRAST(cabins_exp_decay_names) %>% 
  raster::raster() %>% 
  terra::rast()

## Warning in getProjectionRef(x, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = OVERRIDE_PROJ_DATUM_WITH_TOWGS84, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
##  but +towgs84= values preserved

## Warning in rgdal::GDALinfo(fname, silent = ignore.stderr): statistics not
## supported by this driver

## /data/grass/ETRS_33N/u_bb_cuminf/cellhd/private_cabins_sub_inf_nearest_exp_decay1000 has GDAL driver GRASS 
## and has 361 rows and 478 columns

## Warning in getProjectionRef(x, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = OVERRIDE_PROJ_DATUM_WITH_TOWGS84, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
##  but +towgs84= values preserved

terra::plot(c(cabins_exp_decay_r, cabins_exp_decay_grass), 
            main = paste("Exp-decay influence from private cabins - ", 
                         c("R", "GRASS")))

Bartlett decay

Example with ZoI = 1000m

# R
cabins_bartlett_r <- calc_influence_nearest(cabins, transform = "bartlett",
                                            zoi = 1000, where = "R")

# GRASS
name_var <- "private_cabins_sub"
cabins_bartlett_names <- calc_influence_nearest(name_var, transform = "bartlett",
                                                zoi = 1000, where = "GRASS",
                                                quiet = F, overwrite = T)

## projection: 1 (UTM)
## zone:       33
## datum:      etrs89
## ellipsoid:  grs80
## north:      6658900
## south:      6622800
## west:       146900
## east:       194700
## nsres:      100
## ewres:      100
## rows:       361
## cols:       478
## cells:      172558
## Calculating Euclidean distance...
## Reading raster map <private_cabins_sub>...
##    0%   3%   6%   9%  12%  15%  18%  21%  24%  27%  30%  33%  36%  39%  42%  45%  48%  51%  54%  57%  60%  63%  66%  69%  72%  75%  78%  81%  84%  87%  90%  93%  96%  99% 100%
## Writing output raster maps...
##    0%   3%   6%   9%  12%  15%  18%  21%  24%  27%  30%  33%  36%  39%  42%  45%  48%  51%  54%  57%  60%  63%  66%  69%  72%  75%  78%  81%  84%  87%  90%  93%  96%  99% 100%
## [1] "if(A <= 1000.000000, 1 - (1/1000.000000) * A, 0)"
## Calculating Bartlett influence...
## projection: 1 (UTM)
## zone:       33
## datum:      etrs89
## ellipsoid:  grs80
## north:      6658900
## south:      6622800
## west:       146900
## east:       194700
## nsres:      100
## ewres:      100
## rows:       361
## cols:       478
## cells:      172558

# visualize
cabins_bartlett_grass <- readRAST(cabins_bartlett_names) %>% 
  raster::raster() %>% 
  terra::rast()

## Warning in getProjectionRef(x, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = OVERRIDE_PROJ_DATUM_WITH_TOWGS84, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
##  but +towgs84= values preserved

## Warning in rgdal::GDALinfo(fname, silent = ignore.stderr): statistics not
## supported by this driver

## /data/grass/ETRS_33N/u_bb_cuminf/cellhd/private_cabins_sub_inf_nearest_bartlett1000 has GDAL driver GRASS 
## and has 361 rows and 478 columns

## Warning in getProjectionRef(x, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = OVERRIDE_PROJ_DATUM_WITH_TOWGS84, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
##  but +towgs84= values preserved

terra::plot(c(cabins_bartlett_r, cabins_bartlett_grass), 
            main = paste("Linear-decay influence from private cabins - ", 
                         c("R", "GRASS")))

Half-normal (Gaussian) decay influence

Example with ZoI = 1000m

# R
cabins_hnorm_decay_r <- calc_influence_nearest(cabins, transform = "half_norm",
                                               zoi = 1000, where = "R")

# GRASS
name_var <- "private_cabins_sub"
cabins_hnorm_decay_names <- calc_influence_nearest(name_var, transform = "half_norm",
                                                   zoi = 1000, where = "GRASS",
                                                   quiet = T, overwrite = T)

## Calculating Euclidean distance...
## [1] "1.000000 * exp(-0.000003 * pow(A, 2))"
## Calculating half-normal decay influence...

# visualize
cabins_exp_decay_grass <- readRAST(cabins_hnorm_decay_names) %>% 
  raster::raster() %>% 
  terra::rast()

## Warning in getProjectionRef(x, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = OVERRIDE_PROJ_DATUM_WITH_TOWGS84, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
##  but +towgs84= values preserved

## Warning in rgdal::GDALinfo(fname, silent = ignore.stderr): statistics not
## supported by this driver

## /data/grass/ETRS_33N/u_bb_cuminf/cellhd/private_cabins_sub_inf_nearest_half_norm1000 has GDAL driver GRASS 
## and has 361 rows and 478 columns

## Warning in getProjectionRef(x, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = OVERRIDE_PROJ_DATUM_WITH_TOWGS84, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
##  but +towgs84= values preserved

terra::plot(c(cabins_hnorm_decay_r, cabins_exp_decay_grass), 
            main = paste("Halfnormal influence from private cabins - ", 
                         c("R", "GRASS")))

# Euclidean distance
# name_var <- "private_cabins_sub"
# execGRASS("r.resamp.filter", input = name_var, output = "test_neighborhood", filter = "gauss,box", radius = c(500,1500))
# 
# execGRASS("r.mapcalc", expression = "private_cabins_binary = if(isnull(private_cabins_sub), 0, 1)", flags = "overwrite")
# 
# execGRASS("r.neighbors", input = "private_cabins_binary", output = "test_neighborhood", size = 5, flags = "overwrite")
# 
# Zone of Influence
# cumulative zone of influence
# cumulative affected area
# landcape effect, cumulative weight - "Miguet et al"
# Perceived disturbance
# 
# multiply two layers referreing to two different scales of a variable density
# 
# # visualize
# cabins_dens <- readRAST("test_neighborhood") %>% 
#   raster::raster() %>% 
#   terra::rast()
# 
# terra::plot(cabins_dens, main = "Density")
# 
# threshold = 0.05
# execGRASS("r.mapcalc", expression = "areas = if(test_neighborhood > 0.05, 1, 0)", flags = "overwrite")
# 
# cabins_areas_affected <- readRAST("areas") %>% 
#   raster::raster() %>% 
#   terra::rast()
# 
# terra::plot(cabins_areas_affected, main = "areas affected")
# 
# #
# execGRASS("g.extension", extension = "r.area")
# 
# execGRASS("r.clump", input = "areas", output = "patches", flags = "overwrite")
# 
# cabins_areas_patches <- readRAST("patches") %>% 
#   raster::raster() %>% 
#   terra::rast()
# 
# terra::plot(cabins_areas_affected, main = "patch areas affected")